#load data sheet
data <- read.csv("rabbit_data.csv")

#load libraries
library(MCMCglmm)
library(rptR)

#re-code open_retr as binary
data$open_retr[!is.na(data$open_retr)] <- 1
data$open_retr[is.na(data$open_retr)] <- 0

#convert data to long format
ind <- c()
expl <- c()
laten <- c()
month <- c()
open_retr <- c()
labyr_exit <- c()
labyr_retr <- c()
sex <- c()
for(i in 1:nrow(data)){
  ind <- c(ind, rep(as.character(data$ind[i]), 2))
  expl <- c(expl, data$expl_1[i], data$expl_2[i])
  laten <- c(laten, data$laten_1[i], data$laten_2[i])
  month <- c(month, 1, 2)
  open_retr <- c(open_retr, rep(data$open_retr[i], 2))
  labyr_exit <- c(labyr_exit, rep(data$labyr_exit[i], 2))
  labyr_retr <- c(labyr_retr, rep(data$labyr_retr[i], 2))
  sex <- c(sex, rep(data$sex[i], 2))
}

#scale long variables
expl <- scale(expl)
laten <- scale(laten)
labyr_exit <- scale(labyr_exit)
labyr_retr <- scale(labyr_retr)

#convert to data frame
long_data <- data.frame(ind, sex, expl, laten, month, open_retr, labyr_exit, labyr_retr)
long_data$ind <- as.character(long_data$ind)

#scale the original data for sex differences in cognitive traits
data$expl_1 <- scale(data$expl_1)
data$expl_2 <- scale(data$expl_2)
data$laten_1 <- scale(data$laten_1)
data$laten_2 <- scale(data$laten_2)
data$labyr_exit <- scale(data$labyr_exit)
data$labyr_retr <- scale(data$labyr_retr)

#calculate repeatability of expl and laten
expl_rpt <- rptGaussian(expl ~ (1 | ind), data = long_data, grname = "ind", nboot = 10000, npermut = 100000)
laten_rpt <- rptGaussian(laten ~ (1 | ind), data = long_data, grname = "ind", nboot = 10000, npermut = 100000)

#run bivariate analysis with expl and laten as outcome variables, month, sex, open_retr, labyr_exit, and labyr_retr as fixed effects, and ind as random effect, with inverse gamma priors
biv_ig_prior <- list(R = list(V = diag(2), nu = 1.002),
                     G = list(G1 = list(V = diag(2), nu = 1.002)))
biv_for_eff <- MCMCglmm(cbind(expl, laten) ~ trait - 1 + trait:scale(month, scale = FALSE) + trait:scale(sex, scale = FALSE) + trait:scale(open_retr, scale = FALSE) + trait:labyr_exit + trait:labyr_retr,
                        random = ~ us(trait):ind,
                        rcov = ~ us(trait):units,
                        data = long_data,
                        family = c("gaussian", "gaussian"),
                        burnin = 1000, nitt = 100000, thin = 10,
                        prior = biv_ig_prior)

#correlation between expl and laten
expl_laten_corr <- biv_for_eff$VCV[, "traitexpl:traitlaten.ind"]/(sqrt(biv_for_eff$VCV[, "traitexpl:traitexpl.ind"])*sqrt(biv_for_eff$VCV[, "traitlaten:traitlaten.ind"]))
mean(expl_laten_corr)
HPDinterval(expl_laten_corr)

#run univariate analysis with open_retr as outcome variable, sex as fixed effect, and ind as random effect, with inverse gamma priors
univ_ig_prior <- list(R = list(V = 1, nu = 0.002),
                      G = list(G1 = list(V = 1, nu = 0.002)))
univ_sex_open_retr <- MCMCglmm(scale(open_retr, scale = FALSE) ~ scale(sex, scale = FALSE),
                               random = ~ ind,
                               data = data,
                               family = "threshold",
                               burnin = 1000, nitt = 100000, thin = 10,
                               prior = univ_ig_prior,
                               trunc = TRUE)

#run univariate analysis with labyr_exit as outcome variable, sex as fixed effect, and ind as random effect, with inverse gamma priors
univ_sex_labyr_exit <- MCMCglmm(labyr_exit ~ scale(sex, scale = FALSE),
                                random = ~ ind,
                                data = data,
                                family = "gaussian",
                                burnin = 1000, nitt = 100000, thin = 10,
                                prior = univ_ig_prior,
                                trunc = TRUE)

#run univariate analysis with labyr_retr as outcome variable, sex as fixed effect, and ind as random effect, with inverse gamma priors
univ_sex_labyr_retr <- MCMCglmm(labyr_retr ~ scale(sex, scale = FALSE),
                                random = ~ ind,
                                data = data,
                                family = "gaussian",
                                burnin = 1000, nitt = 100000, thin = 10,
                                prior = univ_ig_prior,
                                trunc = TRUE)

#run bivariate analysis with labyr_retr and labyr_exit as outcome variables, sex as fixed effect, and ind as random effect, with inverse gamma priors (to assess correlation between labyr_retr and labyr_exit)
biv_for_cog_corr <- MCMCglmm(cbind(labyr_retr, labyr_exit) ~ trait - 1 + trait:scale(sex),
                         random = ~ us(trait):ind, #G-structure
                         rcov = ~ us(trait):units, #R-structure
                         data = data,
                         family = c("gaussian", "gaussian"),
                         burnin = 1000, nitt = 100000, thin = 10,
                         prior = biv_ig_prior,
                         trunc = TRUE)

#correlation between labyr_retr and labyr_exit
retr_exit_corr <- biv_for_cog_corr$VCV[, "traitlabyr_retr:traitlabyr_exit.ind"]/(sqrt(biv_for_cog_corr$VCV[, "traitlabyr_retr:traitlabyr_retr.ind"])*sqrt(biv_for_cog_corr$VCV[, "traitlabyr_exit:traitlabyr_exit.ind"]))
mean(retr_exit_corr)
HPDinterval(retr_exit_corr)
